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ABSTRACT 

In  this  paper  we  develop  time-domain  state  space  models  for  lossless 
layered  media  which  are  described  by  the  wave  equation  and  boundary  condi- 
tions, Our  models  are  for  non-equal  one-way  travel  times;  hence,  they  are 
more  general  than  existing  models  of  layered  media  which  are  usually  for 
layers  of  equal  one-way  travel  times.  We  develop  state  space  models  for 
two  cases:  (1)  source  and  sensor  at  the  surface,  and  (2)  source  and  sensor 
in  the  first  layer. 
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I.  INTRODUCTION 


In  this  paper  we  develop  time -domain  state  space  models  for  lossless 
layered  media  which  are  described  by  the  wave  equation  and  boundary  condi- 
tions. We  are  specifically  interested  in  models  for  a horizontally  stratified 
nonabsorptive  earth  with  vertically  travelling  plane  compressional  waves, 
and  shall  consider  the  two  cases  depicted  in  Figure  1;  that  is  to  say,  we 
shall  develop  state  space  models  for  K-layer  media  systems  in  which  the 
source  and  sensor  are  located  either  at  the  surface  (the  usual  case, 
considered,  for  example  in  Refs.  1-8)  or  in  the  first  layer  (as  is  the  case 
when  the  first  layer  is  water). 

In  Figure  1 we  characterize  each  layer  by  its  one  way  travel  time, 

T. , velocity,  V\  , and  normal  incidence  reflection  coefficient,  r^ 

(i=l,  2,  . . . , K).  Additionally,  interface-0  denotes  the  surface  and  is 
characterized  by  reflection  coefficient  r^.  Finally,  m(t)  and  y(t)  denote 
the  input  (e.g.,  seismic  source  signature  from  dynamite,  airgun,  etc.) 
to  and  the  output  (e.  g.  , ideal  seismogram  as  measured  by  a geophone  or 
hydrophone)  of  the  K-layer  media  system. 

An  important  use  of  a model  of  a K-layer  media  system  is  to 
generate  synthetic  seismograms;  i.  e.  , to  generate  y(t)  for  a given  m(t). 

This  synthetic  data  can  then  be  used  either  for  preliminary  testing  and 
evaluation  of  a signal  processing  technique  (e.g.  , deconvolution)  or  for 
interpretation  purposes  (Ref.  18).  These  models  are  also  useful  for 
developing  inverse  procedures  by  which  important  parameters,  such  as 
reflection  coefficients  and/or  travel  times,  can  be  extracted  from  measured 
data. 
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Our  time -domain  state  space  models,  as  will  be  seen  below,  are 
quite  different  from  the  more  familiar  z-transform  transfer  function  models 


which  have  appeared  in  the  Geophysics  literature  (Refs.  1-8).  In  the 
Geophysics  literature,  the  assumption  of  equal  one-way  travel  times  is  usually 
made.  Layers  of  different  travel  times  are  built  up  by  inserting  layers  whose 
reflection  coefficients  are  zero.  Our  state  space  models  are  for  non-equal 
one-way  travel  times,  but  can  also  be  applied  to  the  equal  travel  time  case. 

Why  are  we  interested  in  a different  class  of  models  for  what  appears 
to  be  a well  studied  system?  As  is  well  known,  there  is  a vast  literature 
associated  with  systems  which  are  described  by  time-domain  state  space 
models.  Most  recent  results  in  estimation  and  identification  theories,  for 
example,  require  a state  space  model.  These  time-domain  techniques  have 
proven  very  beneficial  outside  of  the  geophysics  field,  and,  we  feel  should 
also  be  beneficial  in  the  geophysics  field.  In  fact,  our  ultimate  objective  is 
to  apply  those  theories  to  the  layered  media  problem;  but,  to  do  so  of  course 
requires  state  space  models. 

One  might  argue  that  it  should  be  possible  to  go  directly  from  the 
z-transform  transfer  function  models,  already  developed,  to  equivalent 
state  space  models.  In  most  cases  this  is  not  practical  since  closed-form 
expressions  for  the  transfer  functions  (e.  g.  , reflection  transfer  function 
Y(z)/M(z))  are  not  available.  Those  transfer  functions  must  be  computed 
from  a set  of  equations  which  are  solved  in  a recursive  manner.  Additionally, 
those  transfer  functions  appear  to  be  limited  by  the  equal  one-way  travel 
time  assumption  (Refs.  2-4  , for  example)  and,  they  appear  to  only  have 
been  published  for  the  case  of  source  and  sensor  at  the  surface  (Figure  la). 


3 


In  this  paper  we  develop  state  space  models  for  both  cases  depicted 
in  Figure  1.  The  case  of  source  and  sensor  at  the  surface  is  treated  in 
Section  II,  whereas  the  case  of  source  and  sensor  in  the  first  layer  is 
treated  in  Section  HI.  Our  state  equations  turn  out  to  be  continuous -time 
equations  with  multiple  time  delays,  and,  are  referred  to  as  causal 
functional  equations.  There  does  not  appear  to  be  any  literature  on  this 
class  of  equations;  hence,  we  also  describe  some  of  their  more  important 
properties  and  two  methods  for  their  computer  simulation  in  Section  II.  A 
connection  between  our  state  space  model  and  transfer  function  models  is 
also  given  in  that  section.  In  Section  III  we  distinguish  between  the  cases 
when  the  sensor  is  either  above  or  below  the  source  and  develop  models 
for  both  cases. 

Some  of  the  material  which  we  discuss  in  Section  II  was  first 
presented  in  Refs.  11,  12,  and  13. 


II.  A STATE  EQUATION  MODEL:  SOURCE  AND  SENSOR 
AT  THE  SURFACE 


A.  State  and  State  Space 

The  starting  point  for  our  developments  is  the  assumption  that  wave 

motion  in  each  layer  is  characterized  by  two  signals  traveling  in  opposite 

directions.  This  assumption  is  a consequence  of  the  lossless  wave  equation. 

th 

Symbols  u^(t)  and  d^(t)  denote  the  upgoing  and  downgoing  waves  in  the  k 
layer,  respectively  (Figure  2).  We  shall  refer  to  u,^(t)  and  d^(t)  as  states. 

Since  the  notions  of  states  and  state  space  may  be  new  to  many 
Geophysics  readers,  we  give  a brief  review  of  them  next.  Our  discussions 
paraphrase  those  in  Refs.  9 and  10. 

The  state  of  a dynamic  system  at  time  t=tg  is  the  amount  of  infor- 
mation at  tg  that,  together  with  the  inputs  defined  for  all  values  of  t>tg, 
determines  uniquely  the  behavior  of  the  system  for  all  t>  t^. 

For  our  layered  media  system,  depicted  in  Figure  1,  the  states 
consist  of  a finite  number  of  variables,  u^(t),  d^(t)  u^ft),  d.^( t),  . . . , 
and  dpi4)-  The  state  of  our  system  can  then  be  represented  by  a column 
vector  x called  the  state  vector,  whose  dimension  is  2Kxl.  The  components 

of  x are  called  state  variables;  hence,  u^(t),  d^(t) u^t),  and  dj^t)  are 

state  variables.  Because  our  input  m(t)  is  real-valued  and  our  state  vector 
is  finite  dimensional  our  state  space,  which  is  defined  as  a 2K-dimensional 
space  in  which  u (t),  d (t),  , u^(t),  d^(t)  are  coordinates,  is  the  familiar 

finite -dimensional  real  vector  space.  The  state  at  time  t of  our  system  will 
be  defined  by  2K  equations  and  can  then  be  represented  by  a point  in 
2K-dimensional  state  space. 


A set  of  state  variables  can  be  associated  with  a given  system  in 


many  ways.  In  other  words  there  exist  a number  of  different  sets  of  state 
variables  for  a given  system  for  a given  input  set.  Which  of  these  is  most 
relevant  depends  on  the  individual  situation,  that  is,  the  nature  of  the 
problem,  the  nature  of  the  input  set,  etc.  We  have  chosen  upgoing  and 
downgoing  signals  in  each  layer  to  be  state  variables  since  these  variables 
are  most  frequently  used  by  geophysicists  and  also  seem  to  be  the  most 
natural  ones  to  use  based  on  wave  equation  theory. 

By  a state  space  model  we  mean  the  set  of  equations  that  describe 
the  unique  relations  between  the  input,  output,  and  state.  It  is  comprised 
of  a state  equation  and  an  output  equation.  The  state  equation  governs  the 
behavior  of  the  state  vector,  x.  For  continuous -time  dynamical  systems  it 
is  usually  a differential  equation,  whereas  for  discrete-time  dynamical 
systems  it  is  usually  a difference  equation.  The  output  equation  relates  the 
output  (or,  outputs)  to  the  state  vector  and  input.  An  example  of  a continuous - 
time  output  equation  is  y(t)  a h'x(t)  + m(t).  In  this  equation  ( )'  denotes 

the  transpose  of  ( ). 

B.  Interface  Equations 

The  starting  point  for  derivation  of  our  state  space  model  is  the 
Figure  3 ray  diagram.  As  in  Refs.  2 and  3,  we  shall  find  it  convenient  to 
draw  ray  diagrams  with  time  displacement  along  the  horizontal  axis,  so  that 
rays  appear  to  be  at  non-normal  incidence  and  so  do  not  overlap  one  another. 
Symbols  u^  and  dj^  denote  the  upgoing  and  downgoing  waves  in  the  layer, 
respectively;  and,  we  adopt  the  convention  that  waves  at  the  top  of  a layer 
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occur  at  present  time  t.  Our  initial  development  is  in  terms  of  u^  and  dj^. 

In  Paragraph  C we  will  find  it  more  convenient  to  work  with  u,  and  d,  . 

k k 

From  Figures  2 and  3 we  see  that  d is  just  a time  delayed  version  of  d' ; 

K KL 

i.  e.  , dk(t)  = 

As  stated  by  Robinson  (Ref.  3),  "the  solution  of  the  wave  equation 

at  each  interface  leads  to  the  definition  of  a reflection  coefficient  r 

J 

associated  with  that  interface.  . . . the  reflection  coefficient  r.,  which  must 

J 

satisfy  | r ^ | < 1,  has  these  properties.  A downgoing  wave  of  amplitude  A in 
layer  j,  upon  striking  interface  j,  is  both  reflected  and  transmitted.  The 
reflected  portion  is  an  upgoing  wave  of  amplitude  r^A  in  layer  j,  so  r. 
represents  the  reflection  coefficient.  The  transmitted  portion  is  a down- 
going wave  of  amplitude  (l  + r.)A  in  layer  j+1,  so  1+r^  represents  the 
transmission  coefficient.  An  upgoing  wave  of  amplitude  B in  layer  j+1  is 
both  reflected  and  transmitted  when  it  strikes  interface  j.  The  reflected 
portion  is  a downgoing  wave  of  amplitude  -r.B  in  layer  j + 1,  and  the 
transmitted  portion  is  an  upgoing  wave  of  amplitude  (l-r^)B.  Hence  -r. 
and  ( 1 - r ^ ) represent,  respectively,  the  reflection  coefficient  and  trans- 
mission coefficient  for  the  upgoing  wave.  These  properties  are  summarized 


in  Table  1 (2).  " 


Table  1.  Reflected  and  Transmitted  Portions 


Reflected 

Portion 

T ransmitted 

Portion 

Downgoing  wave 

Upgoing  wave 

Downgoing  wave 

A 

in  layer  j 

. \sA  . 

in  layer  j 

(l+rj)A 
in  layer  j + 1 

Upgoing  wave 

Downgoing  wave 

Upgoing  wave 

B 

in  layer  j+1 

-rjB 

in  layer  j+1 

(lTj)B 

in  layer  j 

Waveform  u (t+r  ) (Figure  3)  is  made  up  of  two  parts,  namely  the  part  due 

K xC 

to  the  reflected  portion  of  d'(t-T  ) and  the  part  due  to  the  transmitted 

K K 

portion  of  u (t).  It  satisfies  the  equation 

K ri 

Uk(t+V=  rkdk(t-V +(1-rk)uk+i(,)- 

In  a similar  manner,  waveform  d^+j(t)  satisfies  the  equation 

dk+i(t>  ’ '‘W'V  • Vk+i'*'-  <2> 

We  refer  to  Eqs.  (1)  and  (2)  as  the  interface  equations.  These  equations 
are  the  starting  point  for  transfer  function  models,  which  are  very  popular 
in  the  Geophysics  literature  (Ref.  2),  and,  they  are  also  our  starting  point. 

C.  A State  Space  Model 

A state  space  model  for  our  K layer  media  system  is  obtained 

directly  from  Eqs.  (1)  and  (2),  which  are  applicable  at  interfaces  1 through 

K-l  (i.  e.  , for  k=l,  2,  . . . , K-l),  and  comparable  equations  at  the  surface 
th 

and  K interface.  At  the  surface  (Figure  4a),  we  obtain 


y(t)  = rQm(t)  + (l-rQ)  u^t) 

(3) 

d'j(t)  = ( 1 + rQ)  m(t)  - rQ  u^t) 

(4) 
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and,  at  the  K interface  we  assume  that  u (t)  = 0 to  obtain  (Figure  4b) 

KtI 


K<1  + Tk’  = rKdk<t_TKl 

(5) 

= (1  + rK,dk<l-V 

(6) 

Signal  y^t)  in  Eq.  (3)  is  the  measurable  system  output;  hence,  Eq,  (3)  is  the 
outnut  equation.  Signal  d'  (t)  is  also  a system  output;  but,  since  it  cannot 

K.TI 

be  measured,  we  shall  ignore  it  in  following  analyses. 

It  is  convenient  to  group  Eqs.  (1),  (2),  (4),  and  (5)  in  a layer 
ordering,  as  follows. 


d^t)  = -r0u1(t)  + (l  + r0)m(t) 

Ul(t+T  1 ) = rjd'^t-Tj)  + (l-r^u^t) 

dj(«)  =(l  + t._i)d!_1(t-Tjl).r._iU.(t)| 

u.(t+r.)  = r.d!(t-r.)  + (l-r.)u  (t) 

J J J J J J J+l 

dk(t>S(1  + rK-l)dk-l<t-TK-1)-rK-lUK't) 


j=2,  3 K-l 


UK(t  + TK)  = rKdK(t"TK)' 


(7) 


This  system  of  2K  equations  is  not  in  a useful  state  equation  format,  yet, 
since  signals  in  its  left-hand  side  occur  at  t and  delayed  times,  and  signals 
on  the  right-hand  side  occur  at  t,  t-i\  ^ and  t_Tj*  ^ order  to  put  Eq.  (7) 
into  a useful  state  equation  format,  let 

d.(t)  A d!(t-r.)  (8) 


for  all  j=l,  2,  . . . , K.  Observe,  from  Figs.  3 and  2,  that  the  downgoing 
states  d.(t)  occur  at  the  bottom  of  a layer.  Equation  (7)  becomes 
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d l<t+  t 1 ) = - r^^t)  + ( l + rQ)  m(t) 


u^+Tj)  = r1ti1(t)  + (l-rl)  u2(  t) 
dj(t+Tj)  = (l  + r._  . )dj_j(t)  - r.^u.ft) 
Uj(t+Tj)  = r.d.(t)  + ( 1 -r.)u.+1(t) 


j = 2,  3,  . . . , K-  1 


V'^K1  = (ltrK-l)dK-l(,)  ' rK-lUK(t) 


UK,,  + Tk’  = rKdK<‘) 


(9) 


By  means  of  transformation  (8)  each  pair  of  equations  in  (7)  now  only  involves 
two  time  points,  t+T.  and  t.  Equations  (9)  and  (3)  together  represent  a state 
space  model  which  we  refer  to  as  the  layer-ordered  (L-O)  model  in  the 
sequel. 


In  order  to  complete  the  description  of  our  state  space  model  we 
must  specify  initial  condition  information.  The  situation  of  interest  is  the  one 
for  which  all  states  have  zero  values  prior  to  application  of  input  m(t). 
Consider  the  k*^1  layer  (Figure  2).  Then  d^ft)  is  equal  to  zero  until  t - 
Tj  + + • . • +2  T^.  These  facts  are  true  for  all  k=l,  2,  . . . , K;  i.  e.  , 


d.(t)  = 0 V 


r J \ 

t e 0,  E t. 

L i=l  x> 


and 


u.(t)  = 0 v t e 


j 

0,  E T.+  T.  ) 

■ i=l  1 " 


where  j = l,  2,  ....  K. 


(10a) 


(10b) 


As  an  example  consider  a two  layer  system  for  which  Eq.  (9) 
becomes: 
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(11) 


dl(t+Tl)  = ■ roui(t^  + (1  + r0) 
u 1 ( t + T 1 ) = Fjdjft)  + ( 1 - r J ) u2(t) 
d2(t  + T2)  = (l  + r^d^t)  - TjU^t) 

U2(t+T2)=  r2d2(t) 

When  we  solve  this  system  of  equations,  we  must  remember  that  d (t)  = 0 
until  t = t j , Uj(t)  = 0 until  t = 2t  ^ , d^(  t)  = 0 until  t = T ^ +T2-  and  u^(  t)  =0  until 

, = V2V 

D.  Properties  of  the  State  Space  Model 

State  Equation  (9)  is  a dynamical  equation  with  multiple  time  delays. 
It  is  not  a differential  equation,  nor  is  it  a finite -difference  equation.  We 
shall  refer  to  it  as  a causal  functional  equation.  It  is  linear  and  time- 
invariant,  and,  as  is  the  case  with  delay-time  systems,  requires  initial 
value  information  over  initial  intervals  of  time.  We  have  not  been  able  to 
find  any  literature  for  causal  functional  equations. 

Equations  (9)  and  (3)  can  be  expressed  in  more  compact  notation  by 

• • 3$C 

introducing  the  following  2Kx2K  matrix  operator  : 

^ i diagfz^z^z^  z2 ZK' Zk)'  (12) 

where  z.  is  a scalar  operator  used  to  denote  a i\  sec,  time  delay 
( i.  e.  , z.f(t)  = f(t-1\)).  Let 

x(t)  = coMu^t),  dj(t),  u2(t),  d2(t) u^(t),  dK(t))  ; (13) 

* 

This  idea  was  first  suggested  to  us  by  Mr.  Michael  Steinberger,  a former 
graduate  student  in  the  Electrical  Engineering  Department,  at  the 
University  of  Southern  California. 
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then,  Eqs.  (9)  and  (3)  can  be  written,  as 

*x(t)  = Ax(t)  + bm(t)  (14) 

y(t)  = £'x(t)  + rQm(t)  (15) 

where  the  explicit  structures  of  A,  b and  £ can  be  deduced  directly  from  the 
former  equations.  Because  we  do  not  need  this  information  here,  we  do 
not  give  explicit  A,  b,  and  c_  struetures  here  for  Eqs.  (14)  and  115)  Tsee 

Ref.  12.  for  such  structures  for  different  orderings  of  the  states  in  x(t)]. 
From  Eqs.  (14)  and  (15)  we  see  that 


x(t)  = ("^  * -A)  1bm(t)  (16) 

and 

y(t)  = [£'(y  1 -A)'lb  + rQ]m(t)  (17) 

These  equations  provide  us  (conceptually,  at  least)  with  the  solution  of  the 
state  equation  and  with  the  output  as  a function  of  the  input. 

The  reflection  transfer  function  of  the  K layer  media  system  is 
obtained  from  the  Laplace  transform  of  Eq.  (17),  and  is 


Y(s) 

M(s) 


-A)'Vro 


(18) 


where  is  obtained  from  ^ by  setting  z . = e l. 

Equation  (18)  suggests  a straightforward  way  to  compute  y(t)  for  an 
arbitrary  m(t).  First  compute  the  system's  impulse  response,  H(s),  where 
[ recall  that  Y(  s)  = H(  s)  when  M(  s ) = 1 ] , 

H(s)  = c,(^'1-A)'1b  + r0  (19) 

Then,  convolve  h(t)  with  m(t)  to  obtain  yft).  Observe  that  the  right-hand  side 

_ g Q 

of  (19)  is  an  infinite  series  each  of  whose  terms  looks  like  ae  , and,  that 
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J ^ 1 “ S 6 * a 

f {ae  p)  = O' 6( t -g );  hence,  h(t)  is  an  infinite  sequence  of  impulse  functions. 
This  fact  makes  the  convolution  of  h(t)  and  m(t)  quite  easy,  since 

m(t)  * [or6(t-0)]  = arm<t-0)  (20) 

Finally,  we  consider  the  special  case  when  ti=T2=’”’  = Tk^T’  *n 
which  case  ^ = zl,  where  z denotes  the  T sec.  time  delay  and  I is  the 
2K  X 2K  identity  matrix.  In  this  case  Eq.  (14)  can  be  written  as 

x(t+  T ) = Ax(t)+bm(t)  (21) 

Let  us  choose  t = kr  and  assume  that  m(t)  only  has  values  at  t = kT,  in  which 
case  Eq.  (21)  reduces  to  the  following  vector  finite -difference  equation: 

x[(k+l)T]  = Ax(kr)  +bm(kT)  (22) 

In  this  special  case  all  of  the  usual  techniques  (Ref.  10,  for  example) 
associated  with  such  equations  can  be  used  to  solve  for  x(k),  after  which  we 
can  compute  y^kT)  from  the  following  sampled  version  of  Eq.  (15): 

y(kT)  = £'x(kT)  + r^mf  kT ) (23) 

We  do  not  choose  to  follow  this  uniform  travel  time/sampled  data  path, 
because  these  assumptions  seem  too  restrictive. 

E.  Computational  Solutions 

In  this  paragraph  we  briefly  describe  two  computational  methods 
for  practical  computer  solution  of  our  state  space  model  in  Eqs.  (9),  (10)> 
and  (3).  More  detailed  descriptions  of  both  methods  are  given  in  Ref.  13. 

Our  first  method  is  a ray  tracing  technique  in  which  we  define 
mapping  rules  to  track  how  a state  propagates  at  an  interface.  The  rules 
are  obtained  by  observing,  in  Eq.  (9),  what  happens  to  downgoing  or  upgoing 
states  and  are  illustrated  for  downgoing  states  in  Figure  5.  The  complete 
set  of  mapping  rules  are; 
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d.(t) 


Uj(t) 


rQ  m(t) 

y(t) 

j = 0 

rjdi(,) 

-4 

U .( t + T .) 

J J 

. , K 

1 (28) 

(l  + r.)d  (t) 

-♦ 

dj+i(t+V 

j=0... 

. , K-l 

( 1 ’r0)u1(t) 

y(t) 

j=l 

uj-i(t+Tj-i) 

j=2... 

. , K j 

j (29) 

dj(t+T\) 

. , K 

In  order  to  use  Eqs.  (28)  and  (29)  we  create  a state  reference  table. 
Each  element  of  the  states  in  this  table,  called  an  event,  is  characterized  by 
its  time  of  occurrence  and  amplitude.  Our  ray  tracing  procedure  is  based  on 
the  observation  that,  at  an  interface,  an  incoming  signal  branches  into 
reflected  and  transmitted  signals.  Every  time  such  a branching  occurs  we 
have  two  event  points  along  the  time  axis  which  are  stored  in  the  state 
reference  table.  We  search  along  the  time  axis  for  a time  at  which  an  event 
has  occurred.  At  that  time  point  we  map  all  d.  and  u.  states  which  change 
(there  are  only  two  such  states)  via  Eqs.  (28)  and  (29).  As  each  event 
branches  into  two  new  events,  the  table  grows  geometrically.  We  proceed 
along  the  time  axis  looking  up  values  of  d^  and  u^  at  event  points,  until  we 
have  covered  the  domain  of  interest.  To  restrict  the  growth  of  the  table, 
we  use  two  tolerance  parameters,  AMIN  and  5T,  and  collapse  events  that 
occur  at  the  same  time  point  before  we  store  them  back  into  the  table. 
Parameter  AMIN  controls  the  amplitude  of  the  computed  event.  If  that 
amplitude  is  less  than  AMIN,  it  is  set  to  zero.  Parameter  6T  controls  the 
time  separation  of  two  events  below  which  they  are  considered  to  occur  at 
the  same  time. 
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Our  second  method  is  a discretization  technique  in  which  we  convert 


Eq.  (9)  into  standard  finite -difference  equations  by  discretization  of  time. 

Since  the  one-way  travel  times  are  in  general  non-uniform  and  not  a multiple  of  one 

another,  we  have  to  choose  a very  small  time  interval,  A,  for  discretization. 

We  choose  A as  the  greatest  common  factor  of  the  one-way  travel  times,  such 

that  T.=  k.A  (j  = l,  2,  ....  K)  where  k.  is  a positive  integer.  Discretization  is 
J J J 

equivalent  to  dividing  a layer  into  equal  sub-layers  and  inserting  interfaces 

whose  reflection  coefficients  are  zero.  We  illustrate  the  discretization  and 

labeling  of  intermediate  states  for  a layer  where  T . = 4A  in  Figure  6.  The 

facts  that  the  intermediate  downgoing  states  6.(t),  67(t),  and  6,(t)  equal 

dyt+A),  dj(t+2A),  anddj(t+3A),  respectively,  and  the  intermediate  upgoing 

states  n.(t),  u (t),  and  ja  (t)  equal  u.(t+A),  u.(t+2A),  andu.(t+3A),  respec- 
ic  1 j J j 

tively,  are  a direct  consequence  of  the  equations  written  at  interfaces  a,  b, 
and  c.  Those  interface  equations  use  the  fact  that  the  intermediate  reflection 
coefficients  are  zero. 

By  the  discretization  technique  each  layer,  say  the  jth,  expands  to  k. 

layers,  and  each  of  the  k.  layers  is  described  by  two  states:  hence,  the 

J K 

dimension  of  the  state  vector  for  the  discretized  system  is  2^2  k..  The 

i~  1 

easiest  way  to  obtain  the  discrete-time  state  equations  for  the  expanded 


system  is  to  imagine  that  our  earlier  state  equation  derivation,  which  was 

for  the  system  in  Figure  la,  is  for  that  same  system;  but,  now  each  of 

the  K layers  is  further  partitioned  into  k^.  layers,  so  that  Figure  la  is 

now  a system  of  2 £ k.  layers,  each  of  uniform  one-way  travel  time, 
i=l  1 

A.  Of  course,  many  of  the  layers  have  zero  reflection  coefficients  at 
respective  interfaces.  The  discrete-time  counterparts  to  Equations 
(9)  and  (3),  for  the  expanded  system,  are; 
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d (kA  + A) 


and 


u^(kA  + A) 

dj(kA  + A) 

u.(kA  + A) 

J 

dL(k&  + A) 

uL(kA  + A) 


-roVkA)  + ^ + rQ)m(kA) 
r^d^(kA)  + (l-r^u^kA) 


(l  + rj_1)dj_j(kA)-r^_1Uj(kA)  j 

r^dj(kA)  + (l-rJu.+1(kA)  j 

(l  + rL-l)dL-l(kA)-rL-lULfkA) 
rLdL(kA) 


5=2,3, 


. . , L-l 


y(kA)  = rQm(kA)  + (l-r^u^kA) 

where 


1 3 0) 


(3i) 


(32) 


The  non-zero  reflection  coefficients  are  Tq,  r^  , r^  + ^ r^.  All 

other  reflection  coefficients  are  zero. 

Equation  (30)  is  solved  in  an  iterative  manner,  for  k=0,  l,  2,  . . . . 
Usual  matrix  methods  for  obtaining  this  solution  (Ref.  10,  for  example) 
are  terribly  inefficient,  due  to  the  large  value  L will  usually  have  (matrix 
multiplications),  and  the  many  zero  reflection  coefficients  (the  transition 
matrix  for  (30)  is  sparse)  which  are  present.  The  ray-tracing  technique 
applied  to  Eqs.  (30)  and  (31)  is  much  more  efficient.  That  technique  is 
much  simpler  in  the  present  discrete -time  case  because  we  do  not  have 
to  search  along  the  time  axis  for  a time  at  which  an  event  has  occurred. 
Events  occur  every  A units  of  time. 


The  discrete  method  is  simple,  easy  to  implement  and  recursive: 
however,  due  to  the  non-uniform  nature  of  the  delay-terms,  we  have  to  use 
a A small  enough  such  that  it  is  a submultiple  of  all  the  delay  terms.  The 
storage  requirement  is  usually  not  excessive.  However,  we  compute  zero 
for  a lot  of  points  on  the  time  axis  where  no  events  occur.  Thus  the  CPU 
time  may  be  larger  than  for  the  ray-tracing  technique. 
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The  ray-tracing  technique  requires  a slightly  more  complex  program 
to  implement  it.  Its  major  disadvantage  is  the  large  storage  requirement 
for  the  state  reference  table.  Since  this  technique  computes  only  the  actual 
event  points,  it  is  very  efficient  and  fast  for  short  observation  times. 

The  actual  choice  of  a particular  method  depends  on  the  input  data 

set. 


F . Recursive  Reflection  Transfer  Function  Relationship 

We  conclude  this  section  by  making  a connection  between  our  state 
space  model  in  Eqs.  (9)  and  (3)  and  transfer  function  models;  however,  our 
connection  will  be  made  for  systems  with  non-equal  one-way  travel  time. 

For  a K-layer  media  system,  states  u^.(t)  and  d^(t)  have  been  defined 
at  the  top  and  at  the  bottom,  respectively,  of  the  jth  layer  (see  Fig.  2).  Let 
R.(s)  denote  the  transfer  function  between  u.  and  d.  at  the  jth  interface;  i.  e.  , 

<^[u  (t+T  )}  ST  U (S) 

Ri(s)  = ■ * 5^  • (3: 

When  j=0,  we  obtain  the  reflection  transfer  function  between  output  y(t)  and 
source  m(t);  i.  e.  , 


R0(s)  = e 


ST0  U0(3)  A _Y[s) 
Dq(s)  = M(s)  ' 


(34) 


s ince  Tq  = 0. 


We  shall  now  develop  a simple  recursive  relationship  between  R.(s) 


and  R . ,( s).  Consider  our  earlier  state  equations  for  u.  and  d.  , , : 

j+1  n j j+1 
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» 


*■ 


Uj(  t + T j ) = rjdj(t)  + (1~r^)uj  + 1(t) 

(35) 

dj+i<t+Tj+i)  = ,1  + rj,dj,t>  ‘ rjuj+i(t) ' 

(36) 

F rom  the  Laplace  transform  of  Eq.  (35),  -we  find 

R.(s)  = r.  + ( l-r.)U..(s)/D.(s) 

J J J J+1  J 

(37) 

Laplace  transform  Eq.  (36)  and  solve  for  D^(s),  to  show  that 

st 

D,(s)  = [r.U.+1(s)  + e J D ( s )] /( 1 + r) 

(38) 

Substitute  Eq.  (38)  into  Eq.  (37),  and  re-arrange  some  terms  in  the  resulting 


expression  to  see  that 


r.  + z.  , , R.  , ,(s) 

R.(s)=  -J Lj—1 ( j=K- 1 , K-2,  . 


J 


, 1.  0) 


(39) 


1 + rjVlRj+l(s) 


which  is  the  desired  result. 

Equation  (39)  can  be  used  to  compute  the  output  of  a K-layer  media 

system  in  a recursive  manner,  beginning  with  a one  layer  system  ( i.  e.  , one 

layer  on  top  of  a basement  layer)  for  which  we  set  j=K-l.  We  then  iterate 

Eq.  (39)  backwards,  setting  j=K-2,  K-3,  ,..,1,0.  In  order  to  compute 

R (s)  we  need  R ( s ) ; but,  R (s)  can  be  obtained  directly  from  the  very 
K"  I K K 

last  state  equation  in  (9),  u (t+T  ) = r d (t),  as 

K K.  K 

RK(s)=rK.  (40) 

In  the  special,  but  widely  studied  case  of  equal  travel  times,  Eqs.  (39) 
and  (40)  simplify  to 
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J 


.■  +A  (s> 

R.(s)  = -1 , j = K-l,K-2 1,0  (41) 

' 1 + rjZ  V<»> 

RK<Sl  = rK 

- 3T 

where  z = e . Equation  (41  )( or  its  disc rete -time  counterpart,  in  which 
Laplace  transfer  functions  are  replaced  by  z-transform  transfer  functions) 
is  a well-known  result  which  can  be  derived  by  widely  different  methods 
(Refs.  5 and  7,  for  example).  Additionally,  these  recursive  relationships 
occur  (Ref.  14)  in  electric  kernal  functions,  magnetotelluric  input  impedance 
functions,  and  electromagnetic  modified  kernal  functions. 

That  Eq.  (41)  generalizes  to  Eq.  (39)  for  non-equal  travel  times  is 
believed  to  be  a new  result. 
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III.  A STATE  EQUATION  MODEL:  SOURCE  AND  SENSOR  IN  FIRST  LAYER 


In  this  section  we  present  a state  space  model  for  which  source  and 
sensor  are  located  in  the  first  layer.  Such  a model  is  most  important  for 
the  case  of  a marine  environment.  We  shall  consider  the  two  situations 
depicted  in  Figure  7. 

One  approach  for  obtaining  the  state  and  output  equations  is  to  redo 
our  Section  II.  C development  in  which  we  now  include  the  source  and 
sensor  in  the  first  layer.  The  resulting  state  equations  are  somewhat 
different  from  Eq.  (9)  in  that  the  source  m(t)  will  appear  in  the  d^,  u^, 
and  d^  state  equations,  since  m(t)  affects  interfaces  0 and  1. 

A second  approach,  the  one  we  shall  take,  is  one  in  which  we  bring 
m(t)  to  the  surface  so  that  we  can  use  our  Section  II.  C state  equations 
directly.  By  this  artifice,  we  will  reduce  the  derivation  of  the  state 
equations  for  source  in  the  first  layer  to  a problem  which  has  already 
been  solved. 

Observe,  in  Figure  8,  that  there  are  two  "real"  rays  associated 
with  m(t).  Ray  (1)  represents  m(t)  as  it  goes  directly  down  into  the  rest 
of  the  system,  whereas  ray  (2)  represents  m(t)  which  reflects  off  the  surface 
interface  and  then  goes  back  down  into  the  system  (i.e.,  the  ghost).  The 
dashed  ray  is  "imaginary";  i.e.  , it  is  needed  for  construction  purposes 
only,  to  bring  the  component  of  m(t)  along  ray  (1)  up  to  the  surface. 

Let 

me(t)  = IT  + r Q'j  [m(t  +Ta)  " r0m(t'Ta)]  (42) 

At  the  surface  (Figure  9),  we  now  obtain  the  following  state  equation  for  d^: 


20 


(43) 


dl(t  + V = _r0Ul(t*  + (l  + ro^me(t) 

or 

d^t  + Tl*  = -r0ul(t)  + + T a)  " " r a)  (44) 

With  m(t)  brought  to  the  surface,  all  other  state  equations  are  exactly 
the  same  as  in  Eq.  (9);  hence,  the  complete  set  of  state  equations  for  a 
source  in  the  first  layer  is  Eq.  (9)  where  m(t)  is  replaced  by  m (t).  The 
following  initial  conditions,  which  make  use  of  the  fact  that  one  component 
£m(t  + T^jJ  of  the  equivalent  surface  source,  me(t),  occurs  prior  to  time 
zero,  should  be  used: 

dl(t)  = OVtE 
d2(t)  = 0 V t e 

• • • • 

dR(t)  = 0 V t e 
Uj(t)  = 0 V t £ 
u2(t)  = 0 V t e 

• • • 

uK(t)  = 0 V t c 

Now  for  the  sensor  equations,  which  are  a bit  trickier  than  the  state 
equations  because  of  the  "imaginary"  ray  path.  Consider  the  case  where 
the  sensor  is  below  m(t)  (Figure  10a).  In  this  case,  the  rays  which  are 
used  to  bring  m(t)  to  the  surface  do  not  pass  through  the  sensor:  hence, 

yb(t)  = Uj(t  + T d ) + dj(t  + rb)  (46) 

When  the  sensor  is  above  the  source  (Figure  10b)  both  construction 
rays  pass  through  the  sensor,  and 


[°.  VTa  = V 

[°,  T c + T2) 

[°’  TC  +r2  + • * * +tk 
[0,  2Tj  - Ta  = Tc  +rx) 

t) 


[0,  Tc  + 2t7) 


fOT  +t  + • • • + T +2t) 
Lu*  c 2 K-l  K' 


(45a) 


(45b) 
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y = u2(t  + Td)  + d:(t  + rb)  + m(t  - (Ta  - Td))  . m(t  + fT  - r )) 

The  term  m(t  - (t^  - r^))  is  due  to  "real"  ray  (2),  which  really  passes 
through  the  sensor,  but  will  not  be  part  of  the  first  two  state  terms  in 
Eq.  (47),  since  they  are  due  to  inputs  at  the  surface.  When  the  surface 
input  m(t  + Ta)/(1  + Tq)  travels  along  the  "imaginary"  ray  on  its  way  back 
down  into  the  system,  it  registers  the  component  (1  + rA)m(t  + t _ T ,)/(!  + r „ 
at  the  sensor.  Because  this  is  an  imaginary  (j.e.,  non-existent)  component, 
it  must  be  removed  from  the  sersed  signal;  hence,  the  appearance  of 
-m(t  + (Ta  - t d ) ) in  Eq.  (47). 

Figure  11  depicts  m(t)  and  m (t)  for  r.  = 1 and  6 different  values  of 

e 0 

Ta>  ranging  from  Ta  = 1 msec  to  Ta  = 6 msec.  If  the  velocity  in  the  first 

layer  is  5,  000  ft/sec  then  this  range  of  T values  corresponds  to  source 

<1 

depths  ranging  from  50  ft  to  300  ft.  Observe  that  the  most  spike-like  input 
occurs  for  t&  = 1 msec.  As  increases,  mg(t)  becomes  larger  in  amplitude 
and  more  oscillatory,  and  eventually  for  r =5  msec  and  6 msec,  distortion 

2L 

sets  in  the  early  portion  of  me(t). 

The  right-hand  side  of  the  dj  state  equation  in  Eq.  (44)  suggests  that 
our  system  with  its  equivalent  surface  input  is  non-causal:  but,  this  is  not  so. 
We  can  solve  Eq.  (44)  for  d^(t),  as 

dl(t)  ~ "r0Ul(t  " Tl)  ~ rom(t  ‘ Ta  ' Tl)  + m(t  + Ta  * Tl': 
but,  T1-T_  = T , so  that 

1 3L  C 

dl<t)  = _roul(t  ■ Tl)  * rom(t  ■ Ta  ' Tl)  + m(t  ■ V 

We  observe  that  m(t)  occurs  to  the  right  of  zero,  which  confirms  the 
causality  of  our  state  equations. 


i'47) 


) 


(48) 


(49) 
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Our  subsurface  source  and  sensor  model  can  be  simulated  by  either 

of  the  two  methods  described  in  Section  II  - E.  The  mapping  rules  for  the 

ray  tracing  technique  are  similar  to  those  in  Eqs.  (28)  and  (29):  but,  for 

the  present  case  m^ft)  replaces  m(t).  In  the  discretization  technique  A 

must  now  be  chosen  as  the  greatest  common  factor  of  T.,  t, t and, 

and  t^j.  In  the  latter  technique,  we  first  compute  the  system's 

pulse  response,  then  compute  u.(t)  and  d.(t)  via  convolution  between  m (t) 

11  e 

and  the  pulse  response,  and  finally  compute  y (t)  or  y,  (t)  by  delaying  u.(t) 
and  dj(t)  according  to  Eqs.  (46)  and  (47),  respectively. 

The  impulse  response  for  a three-layer  system  is  depicted  in  Figure 
12  for  the  case  of  a fixed  sensor  (t^  = 0.08  sec)  and  varying  source  locations 
(t  = .02,  .04,  and  .06  sec).  The  system  parameters  are:  r_  = l.  00, 
r^  = 0.54,  r^  = 0.34,  r^  = 0.12,  = 0.4,  t = 0.125,  and  = 0.107  sec. 

As  the  source  approaches  the  sensor  the  ghost  reflections  take  longer  to 
reach  the  sensor  (see  spacing  of  first  two  pulses)  and  many  of  the  multiples 
from  the  two  source  components  seem  to  be  getting  closer  to  one  another. 
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IV.  CONCLUSIONS 


f 

We  have  developed  state  space  models  for  lossless  layered  media 
which  are  described  by  the  wave  equation  and  boundary  conditions.  Our 
models  are  for  non-equal  one-way  travel  times,  and  are  therefore  more 
general  than  traditional  transfer  function  models,  which  are  usually  for 
layers  of  equal  one-way  travel  times.  Our  state  space  models  represent  a 
new  class  of  equations,  which  we  call  causal  functional  equations.  These 
equations  are  linear,  time  - invar iant,  continuous -time  equations  with  multiple 
time  delays.  The  impulse  response  of  our  system  is  a sequence  of  unequally 
spaced  impulse  functions. 

We  have  developed  our  state  space  models  for  two  cases:  (1)  source 
and  sensor  at  the  surface  and  (2)  source  and  sensor  in  the  first  layer.  These 
models  can  be  used  either  to  generate  synthetic  seismograms  or  to  develop 
inverse  procedures  for  extracting  reflection  coefficients  or  reflection 
coefficients  and  one-way  travel  times.  Some  recent  work  on  the  extraction 
of  reflection  coefficients  for  source  and  sensor  at  the  surface,  but  from 
noisy  data  is  given  in  Ref.  15.  The  extension  of  the  Ref.  15  results  to  the 
important  case  of  source  and  sensor  in  the  first  layer  is  completed  and  will 
be  reported  on  shortly. 

We  also  wish  to  point  out  that  our  state  space  models  have  led  to 
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some  new  theoretical  results  which  provide  us  with  greater  understanding 
of  the  very  complicated  internal  behavior  of  a layered  media  system.  The 
first  of  these  results  (Ref.  16)  is  a decomposition  of  a seismogram  into  a 
superposition  of  primaries,  secondaries,  tertiaries,  etc.  Each  of  the 
constituents  (e.  g.  , primaries,  secondaries)  is  obtained  from  a state  space 
model  of  dimension  2K.  The  second  of  these  results  (Ref.  17)  uses  this 
decomposition  to  quantify  the  reinforcement  phenomena  that  exist 
within  and  between  the  constituents. 

Finally,  we  reiterate  that  our  ultimate  objective  for  developing  state 
space  models  is  to  apply  time -domain  estimation  and  identification  techniques, 
which  have  proven  to  be  very  beneficial  outside  of  the  geophysics  field,  to  the 
layered  media  problem.  We  are  presently  studying  such  applications. 
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1.  System  of  K-layered  media:  fa)  Source  and  sensor  at 
surface,  and  (b)  source  and  sensor  in  first  layer. 


2.  Definition  of  states  in  layer  k. 


3.  Reflected  and  transmitted  waves  at  interface  k.  From 
Eq.  (8),  d^  (t-rk)  £.dk(t). 

4.  Reflected  and  transmitted  waves  at  (a)  surface  and  (b) 
interface  K. 


5.  Transformation  of  d.(t)  into  neighboring  states  at  fa)  the 
surface,  fb)  jth  interface,  and  fc)  basement  interface. 


6.  Illustration  to  show  insertion  of  intermediate  state  variables 
for  the  jth  layer,  when  k^  = 4. 

7.  Source  (*)  and  sensor  (•)  locations  in  layer  1:  fa)  sensor 
above  source  and  fb)  sensor  below  source. 


8.  Bringing  m(t)  to  the  surface. 


9.  Reflected  and  transmitted  waves  at  the  surface. 


10.  Construction  rays  in  relation  to  fa)  sensor  below  source, 
and  (b)  sensor  above  source. 


11.  m ft)  for  different  values  of  t . 

e a 

12.  Impulse  response  for  a 3-layer  system  for  case  of  a fixed 
sensor  (T.  = 0.08  sec)  and  varying  source  locations:  fa) 

T = 0.02dsec,  (b)  t =0.04  sec,  and  (c)  T =0.06  sec. 
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